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Application  of  a  Nystrom  Method  to 
Solving  3-D  Electromagnetic  Interior 
Scattering  Problems 

INTRODUCTION 

The  impetus  of  our  research  is  electromagnetic  (EM)  dosimetry.  We  are  interested  in  de¬ 
termining  the  distribution  of  electromagnetic  field  strength  inside  a  human  organ  or  even 
the  whole  body  in  presence  of  an  externally  generated  electromagnetic  field.  This  field  can 
come,  for  example,  from  a  cellular  telephone  and  one  is  interested  in  assessing  its  biological 
effects.  Mathematically,  this  is  a  problem  of  an  interior  scattering  problem  where  one  solves 
the  Maxwell’s  Equations  for  an  induced  field  inside  a  finite  body. 

For  our  particular  interest  here,  we  are  interested  in  the  simple  case  in  which  the  external 
source  can  be  modeled  as  a  plane  wave  and  for  frequencies  in  the  microwave  range. 

The  interior  scattering  problem  is  a  classical  problem  and  many  solution  methods  have 
been  proposed.  One  popular  class  of  methods  is  the  finite-difference  time  domain  (FDTD) 
methods  in  which  the  Maxwell’s  Equations  axe  discretized  directly  and  solved  [5, 12].  Another 
class  of  methods,  just  as  popular  as  FDTD,  is  Integral  Equation  (IE)  Methods  [7].  Here 
Maxwell's  Equations  are  cast  into  integral  equations  and  solved  using,  for  example,  the 
Method  of  Moments  (MM)  [4]. 

In  this  work,  we  have  concentrated  on  the  IE  approach.  However,  instead  of  using  the 
Method  of  Moments  (MM),  we  used  a  variant  of  the  Nystrom  Method.  The  Nystrom  method 
is  more  direct  and  does  not  require  approximating  the  sought  field  using  basis  functions  as  in 
MM.  The  merit  of  the  Nystrom  Method  as  a  tool  to  solve  electromagnetic  scattering  problems 
is  rather  controversial.  There  are  some  who  have  reported  success  with  this  method  [2,  13], 
but  others  did  not  [1],  One  of  the  goals  of  our  research  is  to  determine  for  ourselves  the 
merit  of  one  variant  of  this  method.  We  believe  our  approach  is  new. 
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FORMULATION 


In  the  integral  equation  formulation  the  field  E(r)  (we  assume  the  time-harmonic  case  with 
time  factor  e~iut  )  is  written  as  the  sum  of  the  incident  field  E*  and  the  scattered  field  E“: 

E(r)  =  E*(r)  +  E5(r)  (2.1) 

where  the  scattered  field  E8  is  given  by  [9] 

E*(r)  =  (I  +  p  VV-)  jf  g(r,  r')  P(r')  dV  (2.2) 

where 

r(r)  E(r) 
k2(r)~kl 

ejkar 

4ir  r 

|  r  -  r'  | 

The  two  “differentiations”  outside  the  integral  in  equation  (2.2)  are  numerically  undesirable. 
It  would  be  desirable  to  move  the  two  “differentiations”  through  the  integral  and  obtain 

E8(r)  =  [  G(r,  r')  F(r')  dV' 

J  V 

G(r.r')  =  (I+pVV-)j(r,r') 

Unfortunately,  this  is  mathematically  incorrect,  els  the  resulting  integral  is  now  divergent. 
One  can,  however,  still  move  one  “differentiation”  safely  through  the  integral.  In  MM, 
one  can  eliminate  the  remaining  “differentiation”  outside  the  integral  by  moving  it  over 
to  the  testing  functions  [3].  Without  using  the  MM,  we  can  move  it  through  the  integral 
only  indirectly,  resulting  in  a  more  complicated  equation.  In  particular,  by  avoiding  the 
singularity  and  applying  “integration  by  parts”  and  the  divergence  theorem  several  times, 
one  can  arrive  at  the  following  equivalent  formulation  (assuming  F  satisfies  a  local  Holder 
condition): 


F(r)  = 
r(r)  = 

<7(r>  r')  = 


E8(r) 


=  jv  s(r,r')  F(r')  <iV' 

+  h  i  V'V's(r.r')  F(r')  iV 
kl  Jv-v,(r) 
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(2.3) 


+  ^  /,.  VVS(r,r)(F(r')-P(r))il" 

«o  t/Vc(r) 

+  4  /  V'<7(r,r')n'dS'F(r) 

,/ave(r) 

where  K(r)  is  a  sub- volume  of  V  containing  r. 

Two  common  methods  for  evaluating  equation  (2.3)  exist,  depending  on  the  choice  of 
Ve  (r).  By  decomposing  the  body  V  into  finite  volume  elements  (e.g.  tetrahedra),  one  can 
naturally  take  Ve (r)  to  be  the  volume  element  that  contains  r  [14]. 

The  other  common  method  allows  V^(r)  to  vanish  around  r.  This  leads  to  an  appar¬ 
ently  simpler  integral  equation  [14].  However,  it  requires  the  evaluation  of  a  principal- value 
integral,  which  is  numerically  difficult. 

The  approach  we  have  taken  here  is  to  let  (r)  =  V.  This  leads  to  the  equation 

Es(r)  =  /  <7(r,r')  F(r')  dV'  +  1  /  V'P(r,r')n'  dS'  F(r) 

JV  «0  J QV 

+  J  V'VV(r.r')  (F(r')  -  F(r))  <JV' 

=  A(r)F(r)+/  G(r,  r')(F(r')  -  F(r))  dV1  (2.4) 

JV 

where 

A=4  /  V'g(r,r')n' dS'  +  /  g(r,r')ldV' 

K0  J  av  Jv 


For  a  specified  r,  A  is  completely  known,  so  it  can  be  accurately  estimated. 


THE  NYSTROM  METHOD 

The  Ny strom  Method  assumes  the  integrals  in  Equations  (2.3  or  2.4)  can  be  approximated 
well  by  quadrature  formulas  using  N  quadrature  points  {r;}^ .  If  Equation  (2.1)  is  also 
evaluated  at  these  same  N  points,  then  one  obtains  a  linear  system  of  N  vector  equations 
in  the  N  unknown  vectors  {E(r;)}^:1.  Inverting  this  system  will  yield  the  values  of  {E(r,)}, 
which  in  turn  can  be  used  to  estimate  the  value  of  the  electromagnetic  field  E(r)  at  any 
other  point  r  by  using  a  discretized  version  of  equation  (2.1). 

The  main  sources  of  errors  in  the  Nystrom  method  as  formulated  here  are  (1)  the  sin¬ 
gularity  of  some  of  the  integrals,  (2)  the  numerical  approximation  of  the  integrals,  and  (3) 
the  possible  ill-conditioning  of  the  system  of  equations.  We  will  only  address  the  first  two 
sources  of  errors  in  the  following. 
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In  some  of  our  previous  numerical  experiments  [8]  with  equation  (2.3)  using  MM,  we  have 
encountered  some  numerical  difficulties.  Unsatisfactory  results  were  observed  whenever  the 
contrast  r  is  large  or  when  more  tetrahedra  were  used  to  represent  the  volume  V.  This  could 
be  caused  by  a  multitude  of  things  including  the  improper  handling  of  the  singularity,  the 
poor  choice  of  basis  functions  and/or  testing  functions,  inappropriate  dot-product,  low-order 
numerical  integration,  sub-optimal  choice  of  tetrahedra,  or  poor  meshing  in  general. 

By  using  the  Nystrom  method,  we  eliminate  the  need  to  consider  basis/testing  functions 
and  to  some  extent  (as  we  will  we  see  later)  meshing. 


SINGULARITY:  A  1-D  TEST 

Equation  (2.4)  is  theoretically  correct  if  we  assume  F  satisfies  a  local  Holder  condition  so 
that  the  integral 

Jv  V'V'<7( r,r')  (F(r')  -  F(r))  dV'  (3.5) 

is  convergent,  despite  the  singularity  at  r.  Since  F  is  unknown,  one  way  to  incorporate  this 
assumption  on  F  (and  we  must)  into  a  numerical  quadrature  is  to  assume  the  contribution  of 
the  singular  point  r  to  the  volume  integral  in  Equation  (3.5)  is  negligible.  More  precisely,  the 
volume  integral  over  a  sufficiently  small  volume  surrounding  r  is  negligible.  This  is  basically 
a  reiteration  of  the  fact  that  Equation  (3.5)  is  convergent  when  F  satisfies  a  local  Holder 
condition. 

To  get  an  idea  of  the  consequence  of  this  assumption,  we  studied  the  following  1-D 
problem  analogous  to  the  actual  3-D  problem  described  in  Equation  (2.4). 


where  typically, 


PH*)  -  (j 

+  “  £5  jf  «(*.  *')  /(*')  *')  =  «'"» 


f(x)  =  t{x)u{x),  t(x)  =  k2(x)  —  kl 

0(x)  =  — ,  a=i 


r(xY 


kl 


Furthermore,  we  assume  (again  motivated  by  the  actual  3-D  problem) 

g{x,x')  =  g{x,x') 

9x{x,x)  =  - gx-(x,x' ) 


(3.6) 
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J  9{X:  ^  )  dx 

is  convergent 

1  dx 

exists 

x  —  t 

f  9x,x,{xi  x 
f0 

is  divergent 

Analogous  to  the  3D  case,  this  equation  can  be  reformulated  as 


[P  -  r(x)  ~  J0  9{x>x')dx'}  /(*)  ~  Jq  [g{x,x') 

+  agx.xi(x,x')](f(x')  -  f(x))dx' 

=  ul(x) 

where 

r(x)  -  a\gx(x,  0)  -  gx(x,  1)] 


(3.7) 


(3.8) 


A  function  g(x,x')  that  has  all  the  properties  listed  in  (3.7)  is 

g(x,  x)  =|  x  —  x  |  ( In  |  x  —  x  |  —1 )  (3-9) 

This  is  the  function  we  have  used. 

In  the  Nystrom  Method  we  used,  the  N  quadrature  points  (and  corresponding  weights) 
chosen  to  evaluate  the  integrals  in  the  IE  axe  the  those  determined  by  the  Gauss-Legendre 
quadrature  points  (GLQPs)  [10]  on  [0,1].  More  sophisticated  choice  of  quadrature  points 
(and  weights)  that  depends  on  the  singularity  in  question  could  have  been  designed  [6],  but 
was  not  done  in  this  exploratory  study. 

The  sum  of  the  weights  is  always  equal  to  the  length  of  the  integration  interval.  We  might 
view  the  GLQPs  as  a  partition  of  [0,1]  into  sub-intervals  each  of  which  is  centered  around 
a  GLQP  and  whose  length  is  equal  to  the  corresponding  weight.  The  assumption  that  the 
solution  satisfies  a  local  Holder  condition  implies  that  the  contribution  from  (2.4)  is  finite 
and  negligible  if  the  norm  of  the  partition  is  small.  In  short,  we  may  make  the  following 
approximation: 


Jq  9xx{xi,x')(f(x')  -  f(xi))dx' 

=  £  wjgxx{xi,Xj)  (, f(xj )  -  f{xi))  (3.10) 

j¥* 
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We  generated  simulated  data  by  using  the  function 

f(x)  =  x(l  -  X ) 


(3.11) 


to  calculate  an  “incident  wave”. 

The  results  are  shown  in  Figure  1.  It  shows  the  actual  solution  and  the  calculated  solution 
for  IV  =  10  and  N  =  25. 


Figure  1 .  Calculated  vs  Actual  Field  f(x) 


x 

For  this  example,  the  calculated  solution  approximates  well  the  actual  solution  even  when 
N  is  as  low  as  25.  The  “relative  permittivity”  here  is  3. 
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AN  ITERATED  INTEGRAL  METHOD 


The  apparent  success  of  the  1-D  analog  problem  holds  out  hope  that  the  same  might  be  true 
for  the  actual  3-D  problem.  Since  the  formulation  of  the  problem  in  the  1-D  case  is  identical 
to  that  in  the  3-D  case,  it  remains  only  to  choose  a  3-D  integration  method  to  complete  the 
transition  from  1-D  to  3-D. 

The  following  is  a  brief  summary  of  our  experimenting  with  different  possible  3-D  inte¬ 
gration  methods,  cumulating  in  the  Iterated  Integral  Method  (IIM). 

A.  Integration  on  Tetrahedra 

Methods  of  integration  over  tetrahedra  abound  in  the  literature.  (Notably  the  classic  work 
of  Stroud.)  We  have  implemented  some  of  these  methods  by  approximating  the  body  in 
question  by  a  family  of  tetrahedra  or  a  family  of  curvilinear  tetrahedra  (for  a  better  approx¬ 
imation)  and  using  different  integration  formulas  (differing  in  orders)  on  each  tetrahedron 
or  curvilinear  tetrahedron  (after  appropriate  transformation.)  The  results  are  mixed  but 
generally  disappointing.  We  remark  in  passing  that  by  increasing  the  order  of  the  method 
of  integration  on  each  tetrahedron,  one  invariably  increases  the  number  of  integration  points 
in  each  tetrahedron  and  thereby  increases  the  number  of  unknowns  and  hence  the  size  of 
the  system  matrix.  The  effect  of  increasing  the  order  of  the  method  of  integration  on  the 
accuracy  of  the  final  solution  is  not  clear  at  this  time.  Our  preliminary  studies  also  seem 
to  indicate  an  increase  in  condition  numbers  of  the  system  matrices  as  the  order  of  inte¬ 
gration  is  increased.  This  suggests  a  subtle  difference  between  numerical  integration  of  a 
known  function  and  that  of  an  unknown  function  (as  in  solving  an  IE):  merely  using  a  more 
accurate  numerical  integration  method  does  not  necessarily  yield  a  more  accurate  IE  solver. 


B.  Extrapolation 

One  way  to  increase  the  order  of  the  integration  method  without  increasing  the  number  of 
integration  points  is  to  use  points  outside  the  integration  domain.  So,  instead  of  using  only 
points  of  integration  within  a  tetrahedron,  we  can  also  use  a  pre-specified  number  of  points 
in  neighboring  tetrahedra.  We  have  implemented  one  such  method.  In  this  method,  the 
set  of  integration  points  will  consist  of  the  ’’centers”  of  all  the  tetrahedra  in  the  system. 
(This  is  different  from  the  1-D  case  where  the  integration  points  are  not  predefined,  but 
are  chosen  to  optimize  the  order  of  the  integration  method.)  For  each  tetrahedron,  T,  we 
picked  a  pre-specified  number  of  nearest  (independent)  centers  as  the  integration  points 
for  this  tetrahedron,  i.e.,  weights  associated  with  these  selected  integration  points  and  for 
this  tetrahedron  will  be  calculated  so  that  certain  moments  are  integrated  exactly  on  the 
tetrahedron  T. 
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We  observed  that  the  weights  generated  by  this  method  are  often  negative  and  are  sus¬ 
ceptible  to  round  off  errors  as  well  as  errors  incurred  by  approximating  the  given  body  by 
tetrahedra. 

We  did  not  observe  a  significant  improvement  over  the  previous  methods  in  which  only 
interior  quadrature  points  are  used  in  the  quadrature  formulas  for  approximating  integrals 
over  tetrahedra. 


C.  Iterated  Integral  Method  (IIM) 

Under  the  assumption  that  F  satisfies  a  local  Holder  condition,  all  integrals  in  Equation  2.4 
are  absolutely  convergent  and  we  can  evaluate  them  as  iterated  integrals  (Fubini’s  Theorem.) 
Consequently,  since  each  iterated  (volume)  integral  is  made  up  of  three  separate  1-D  integrals, 
one  cam  easily  generate,  for  example,  3-D  Gauss-Legendre  quadrature  points  and  weights 
from  1-D  Gauss-Legendre  quadrature  points  and  weights.  As  long  as  the  material  boundaries 
(internal  and  external)  of  the  body  in  question  are  known,  the  Gaussian  points  and  weights 
can  be  calculated  rather  straightforwardly.  In  particular,  no  explicit  discretization  (meshing) 
of  the  body  is  required.  This  can  be  a  significant  advantage  this  method  has  over  the  more 
conventional  methods  in  which  meshing  is  required,  as  meshing  is  sometimes  the  Achilles’ 
heels  of  these  methods. 


A  3-D  TEST  OF  NYSTROM  METHOD  USING  IIM 

We  tested  the  Nystrom  Method  using  IIM  (NM-IIM)  on  an  integral  equation  identical  in 
form  to  equation  (2.1)  wherein  equation  (2.4)  is  used  for  the  scattered  field.  However,  the 
test  problem  was  chosen  to  have  no  singularity,  sis  we  were  only  interested,  at  the  initial 
stage,  in  testing  the  NM-IIM  in  a  setting  insulated  from  other  factors  that  may  complicate 
the  problem. 

The  irradiated  body  is  a  5  cm  sphere.  For  simplicity,  we  used 


E(r)  =  [x,y,z] 

9{ r,r')  =  (a :  -  x')A  +  (y  -  y')A  +  (z  -  z')* 

in  the  equation 


B  F(r)  -  /  G(r,  r')  (F(r')  -  F(r))  dV' 

Jy 

-  Einc(r) 
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to  generate  an  “incident  wave”.  Here 

B  =  o( r)I-^  /  V,^(r,r,)n'  dS' 

%  J  qv 

“W  =  ^)-Iv^dV‘ 

P(r)  =  r(r)E(r) 

Figure  2  compares  the  real  part  of  the  z-component  of  the  calculated  field  with  that  of 
the  actual  field  along  the  z-axis.  For  this  simple  test,  this  is  just  z.  Thus,  the  straight  line 
in  Figure  2  is  expected. 
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A  total  of  151  quadrature  points  in  the  sphere  were  picked  (by  the  solver).  The  “relative 
permittivity”  was  chosen  to  be  52.  Judging  from  this  test,  albeit  simple,  the  IIM  seems  to 
be  working. 


Figure  2.  Actual  vs  Calculated  Re(Ez)  (er=52:  N=151) 
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APPLICATIONS 


Finally,  we  applied  the  3D  NM-IIM  to  a  realistic  problem  consisting  of  a  dielectric  sphere 
incident  by  a  1-GHz  plane  wave  propagating  along  the  z-axis  towards  positive  infinity.  The 
radius  of  the  sphere  is  5  (cm)  and  the  conductivity  is  set  to  0.015625  (mhos/m). 

A.  Homogeneous  Sphere 

The  results  of  applying  NM-IIM  to  the  above  problem  using  differing  relative  permittivities 
(er)  axe  shown  in  the  Figure  3.  It  compares  the  reed  part  of  actual  x-component  of  the  E  field 
along  the  z-axis  to  that  calculated  by  NM-IIM.  The  dotted  lines  in  the  Figure  3  represent 
the  actual  (Mie)  solutions,  one  for  each  relative  permittivity.  Here  eT  «  1+,  2+,  3+,  and  13+. 

There  is  generally  a  fairly  good  agreement  between  the  actual  and  the  calculated,  al¬ 
though  again  performance  noticeably  deteriorates  with  increased  relative  permittivity  or, 
equivalently,  with  higher  contrasts.  The  number  of  quadrature  points  used  for  this  sphere  is 
608. 

B.  An  Inhomogeneous  Sphere 

By  manipulating  the  limits  of  integration  in  IIM,  the  NM-IIM  can  also  be  applied  to  inhomo¬ 
geneity  bodies.  More  generally,  one  can  vary  the  grid  size  (resolution)  from  one  sub-region 
to  another  in  the  body  quite  straightforwardly. 

The  Figure  4  shows  the  result  of  using  NM-IIM  on  the  above  sphere  using  different 
resolutions  in  two  different  sub-regions.  The  sphere  is  divided  into  two  sub-regions,  an  inner 
sphere  and  an  outer  layer.  The  radius  of  the  inner  sphere  is  2.5  cm,  whereas  the  complete 
sphere  again  has  radius  5  cm.  The  outer  layer  is  chosen  to  have  twice  the  resolution  as 
that  of  the  inner  sphere.  We  also  compare  this  variable  grid  case  to  a  uniform  grid  case. 
The  resolution  of  the  uniform  grid  is  equal  to  that  of  the  inner  sphere  of  the  variable  grid 
case.  A  total  of  151  points  is  used  for  the  uniform  grid  case  and  612  points  for  the  variable 
grid  case.  The  relative  permittivity  is  approximately  3.  Figure  4  shows  the  calculated  and 
actual  real  part  of  the  x-compoenent  of  the  E  field  along  the  z-axis  of  the  sphere.  Again,  the 
incident  plane  wave,  linearly  polarized  in  the  x-direction,  is  propagating  from  the  negative 
z-direction.  linearly  polarized  in  the  x-direction. 

Again,  the  calculated  field  compares  quite  well  with  the  theoretical  solution  and  the 
added  resolution  in  the  outer  layer  has  noticeably  improved  the  accuracy  of  the  method. 
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Figure  4.  Calculated  Re(Ex)  Using  Variable  Grid  and  Uniform  Grid  vs  Actual  (er  =  3+) 
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DISCUSSION  AND  CONCLUSION 


The  Nystrom  Method  in  combination  with  the  Iterated  Integral  Method  (NM-IIM)  as  imple¬ 
mented  here  has  out-performed  other  methods  that  we  have  previously  tried.  Because  it  is 
“grid-free”  and  relatively  easy  to  implement,  it  can  easily  be  used  on  other  geometry  (besides 
spheres)  like  layered  ellipsoids,  finite  layered  cylinders,  or  even  simple  human  organs. 

There  is  obviously  much  room  for  improvement.  The  first  approximation  we  used  in 
treating  the  singularity  can  obviously  be  improved.  The  fact  that  the  estimates  near  the 
boundary  is  generally  poorer  than  those  in  the  middle  may  be  consequence  of  this  first 
approximation.  In  the  neighborhood  of  a  singularity,  the  integrand  is  a  function  of  the 
Jacobian  of  the  unknown  F(r).  The  assumption  that  this  contribution  is  small  may  not  be 
appropriate  near  a  boundary. 

We  should  probably  use  even  more  sophisticated  1-D  Gaussian  quadrature  formulas  that 
will  take  better  account  of  the  singularity  in  the  problem,  similar  in  flavor  to  that  discussed 
in  [6].  Hopefully,  this  will  then  allow  us  to  address  problems  with  even  higher  contrasts  that 
one  routinely  encounters  in  human  organs. 

A  rigorous  numerical  analysis  of  NM-IIM  is  obviously  needed.  Good  performance  of  a 
method  on  test  problems  only  lends  credence  to  the  method,  but  never  an  absolute  proof. 
The  same  analysis  can  invariably  lead  to  improvement  of  the  method. 

Finally,  as  is  always  the  case,  we  would  like  to  be  able  to  solve  more  practical  and  more 
realistic  problems  (bigger  geometry,  higher  frequencies,  and  severer  inhomogeneity)  and  to 
do  so  faster.  Parallel  computing  on  high  performance  computers  is  the  most  obvious  answer. 
We  have  just  started  to  explore  this  possibility. 
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